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ABSTRACT 



Gravitational instability plays an important role in driving gas accretion in massive pro- 
tostellar discs. Particularly strong is the global gravitational instability, which arises when the 
disc mass is of order 0. 1 of the mass of the central star and has a characteristic spatial scale 
much greater than the disc's vertical scale-height. In this paper we use three-dimensional nu- 
merical hydrodynamics to study the development of gravitational instabilities in a disc which 
is embedded in a dense, gaseous envelope. We find that global gravitational instabilities are 
the dominant mode of angular momentum transport in the disc with infall, in contrast to oth- 
erwise identical isolated discs. The accretion torques created by low-order, global modes of 
the gravitational instability in a disc subject to infall are larger by a factor of several than an 
isolated disc of the same mass. We show that this global gravitational instability is driven by 
the strong vertical shear at the interface between the disc and the envelope, and suggest that 
this process may be an important means of driving accretion on to young stars. 
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1 INTRODUCTION 

Accretion discs play a fundamental role in many aspects of as- 
trophysics. Objects as diverse as planets, stars and super-massive 
black holes are all thought to acquire significant fractions of their 
mass through disc accretion, and consequently understanding ac- 
cretion disc physics is important in understanding the formation 
of all of these objects. Critical to our understanding of accretion 
discs is the process of angular momentum transport, but despite 
many years of research on this subject we still do not fully un- 
derstand the mechanism(s) by which angular momentum is trans- 
ported in gaseous discs. In many cases we believe that magneto- 
hydrodyn amic instabilities, such as the m agnetorotational instabil- 
ity (MRI. lBalbus & HawIevlll99ll , [T99^ is the dominant transport 
mechanism. However, some systems, notably protostellar discs, 
are insufficient ly ionized for the MRI to operate everywhere (e.g., 
lGammiellT99^ . and it is also not clear whether or not the MRI 
can drive accretion rates as high as those which are observed 
jKing. Pringle & Livio 2007). Consequently, it is still desirable to 
investigate other mechanisms for angular momentum transport in 
discs. 

One such mechanism which has received considerable inter- 
est in recent years is angul ar momentum transport by gravitat ional 
instabilities (GIs; see, e.g.. lDurisen et al. l l200ilLodatoll2008l and 
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references therein). Gaseous di scs in Keplerian rotation become un- 
stable to self-gravity when the' Toonir^ jl964) Q parameter is less 
than some critical value of order unity. The Toomre parameter is 
defined as 



Q 



nGL' 



(1) 



where Cj is the sound speed of the gas, Q. is the angular frequency 
and S is the surface density. Shearing discs generally become unsta- 
ble to non-axisymmetric perturbations before axisymmetric ones, 
so GIs in discs initially manifest themselves as spiral density waves. 
It was recognised long ago that such spiral density waves c an trans- 
port angular momentum ij Lvnden-Bell & KalnaisI Il972h but de- 
tailed study of the non-linear development of GIs in gaseous discs 
has only recently become possible. This process has been studied 
in great detail using numerical hydrodynamics, and we now have a 
well-established picture whereby angular momentum transport by 
GIs is primarily governed by disc thermodynamics. GIs in isolated 
thin gas discs tend to evolve to a self-regulating state, where the 



ing (e.g., Gammie 200ll; 


Lodato & Ricell2004 iMeii'a et al.l 


2OO5I; 


ICossins. Lodato & Clarke 


I2OO9I), and althoush eravitv is a 


long- 



range force, "global" effects generally do not dominate unless the 
disc mass is an appreciable fraction (> 25%) of the mass of the cen- 
tral o bject jLaughlin. Korchagin. & Adamsl [19981 : iLodato & Ricd 
I2OO5I) . In this picture the efficiency of angul ar momentum trans- 
port can be parametrized in terms of a classical fshakura & SunvaevI 
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Orbital period / yr 

Figure 1. Maximum sustainable accretion rate in a critically self-grav itating 
disc with a = 0.1, computed following tiie procedure described in Levin' 



2007 ) and using the opacities K{p, T) of Bell & Lin ( 1994) and Bell et al. 



19971) . The sharp jump in the critical accretion rate at an orbital period of 



^ 300yr is caused by the transitio n between th e optically thick inner disc 
and optically thin outer disc (Matz ner & Levinl20 05V This corresponds to a 
radius of =: 40 AU for a IMg central star, or ^ lOOAU for a IOMq star. The 
maximum sustainable accretion rate at larger radii is small, ~ lO^^Moyr"' ; 
in the "local limit", larger accretion rates lead to fragmentation. In some 
cases external irradiation can be the dominant source of heating, impos- 
ing a temperature "floor" (denoted by Tmin) and enhancing the maxi mum 
accretion rate. Similar figures can be found in .Clarke. t2009 ) and Rafikovl 
i2009l) . 



lll973h a-prescription jGammi j200ll : lLodato & Ricel2004) . where 



1 



9r(7- D/cooi" 



(2) 



Here /cool is the local cooling time-scale and y is the adiabatic index 
of the gas. Faster cooling leads to deeper spiral density waves (i.e., 
with higher density contrasts), and thus to more efficient transport 
of angular momentum. However, if the cooling becomes too rapid 
the disc is unable to maintain its self-re gulating state, and the GIs 
instea d lead to fragmentation of the disc ( lGaminiel200ll : lRice et al.l 
|200^ . This in turn imposes a maximum efficiency at which angular 
momentum can be transported by GIs without leading to disc frag- 
menting, and numerical simulations place typically this "fragmen- 
tation boundary" at q-qi £ 0.1 (corresponding to ,8 = fcooiQ ^ 3- 
5, e.g.. lGammi3l200ll : lRic"e et al.ll2003l : lRice, Lodato, & Atmitag^ 



I2OO5I) . When extended to consider discs with realistic opacities, 
these results imply a maximum accretion rate tha t can be sustained 
by GIs in a self-regulating state (e.g., Levin 2003l: lMatzner & LevinI 
l2005ljLevinl200l1ciard2009l : lRafikov.2009. . see also Fig.[D- Ex- 
cept at very small radii this rate is low, ~ 10"*Moyr"', and this 
raises questions as to how many astrophysical objects are able to 
accrete their mass in a plausible time-scale. The star may continue 
accreting bound clumps of gas even after the disc fragments (eg., 
IVorobvov & Basull201(]h . but the details of this process remain un- 
certain. 

To date most numerical studies of GIs have looked at isolated 



self-gravitating discs, but in reality it seems likely that most gravi- 
tationally unstable discs will still be subject to some level of infall 
on to the disc. Indeed, in many cases it is likely that the instanta- 
neous infall rate on to the disc exceeds the accretion rate through 
the disc. For example, observed accretion rates on to protostellar 
discs are typically an order of magnitude la rger than the accre- 
tion rates on to the protostars the mselves (e.g.. lKenvon et aljl990l ; 
ICalvet, Hartmarm & Stroml '2000). Similar discrepancies between 
infall and disc accretion r ates have been f ound in models of low- 
mass star formation (e.g.. lVorobvovll2009l) , and in models of star 
forma tion in black hole accretion discs (e.g.. lMilosavlievic & Loebl 
l2004h . In this paper we present an initial investigation of this 
problem, by using three-dimensional numerical hydrodynamics to 
follow the evolution of a self-gravitating accretion disc subject 
to quasi-spherical infall. In Section |2] we present our numerical 
method, and in Section |3] we discuss the results of our simula- 
tions. We find that infall on to the disc can substantially enhance 
the efficiency of angular momentum transport, through the excita- 
tion of low-order, global, spiral density waves. We discuss the con- 
sequences of this result for real astrophysical systems, along with 
the limitations of our analysis, in Section |4] and summarize our 
conclusions in Section|5] 



2 NUMERICAL METHOD 

Our simulations are conducted using the publicly- available 
smoothed-particle hydrodynamic (SPH) code gadget-2 dSpringell 
l2005h . We have modified the code to include a simple scale-free 
cooling prescription, as used in previous sim ulations jGammiel 
I2OOII : iLodato & Ricel2004l : ICossins et al.l2009l) . which has the fol- 
lowing form: 



dui 
dt 



^cool 



(?) 



Here is the internal energy of particle i, and the cooling time- 
scale tcooi is proportional to the local dynamical time-scale thus 



Operationally, the cooling time-scale is computed as 



tcool - P 



I 



(4) 



(5) 



where Rj is the cylindrical radius of the rth particle. The cool- 
ing time thus depends only on radius, and does not vary with z 
or with the instantaneous orbital speed (which can be perturbed 
significantly in unstable discs). As mentioned in Section[T] previ- 
ous simulations of self-gravitating discs have found that values of 
j8 < 3-5 result in fragmentation of the di sc, while larger values lead 
to transport of angular momentum (e.g.. lGanimiell20()ll ; lRice et al.l 
I2003ll2005h . We do not wish to see disc fragmentation due to rapid 
cooling alone, and therefore set /? = 7.5 throughout. We adopt an 
adiabatic equation of state, with adiabatic index 7 = 5/3. 

We make use of a single sink particle as the central gravitat- 
ing mass, w hich accretes all ga s particles within its sink radius (as 
described in ICuadra et alj2006l) . This is primarily a numerical con- 
venience, used in order to prevent the time-step being limited by 
a small number of SPH particles at very small radii, and has no 
physical effect on the simulations. We use the standard Barnes- 
Hut formalism to calculate the gravitational force tree, and use 
N„gfc = 64±2 as the number of SPH neighbours. We allow a variable 
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gravitational softening length, whic h is equal to th e SPH smoothing 
length throughout (as demanded bv lNelsoi]|2006l) . The simulations 
are scale-free: we use a system of units where the central gravitat- 
ing mass has an initial mass M^, = 1, the inner edge of the disc is 
at /? = and the time unit is the orbital period at R = 1. (Thus 
G = 4n^ in code units.) 



2.1 Artificial viscosity 

We adopt the sta ndard Monaghan-Gingold-Balsara form for the ar- 
tificial viscosity iMonaghan & Gingold l983VB alsaral 19951) . as de- 
scribed in Equations 11-12 of Springel (2005j). This prescription 
contains both linear and quadratic terms (characterised by the pa- 
rameters Q-j,,/, and /3sph respectively, with flsph = 2aspi,), and the 
"Balsara switch" which acts to limit the artificial viscosity in pure 
shear flows. We adopt a,,,;, = 0.3 throughout. 

As we are primarily interested in how angular momentum 
is transported in our simulations, great care must be taken to en- 
sure that this transport is not dominated by numerical effects. It is 
well-known that SPH artificial viscosity can drive significant an- 
gular momentum tran sport in disc simulations (e.g.. iMurravll 1 99^ : 
lLodato&Ric3l2004l) . so we have conducted tests to ensure that 
this is not the dominant source of transport in our models. From 
our standard disc initial conditions (see Section |2.2. 1 1 below) we 
ran simulations with the self-gravity of the gas turned off; with this 
set-up, Reynolds stresses due to numerical effects (primarily the ar- 
tificial viscosity) are the only source of angular momentum trans- 
port. By_£X2ressmgthis stress in units of the local pressure (see, 
e.g.. lLodato & Ricel2004) we can parametrize the efficiency of the 
numerical transport as a familiar a-parameter thus: 



2 dv,-Sv^ 

3^r 



(6) 



where S\ = v- (v) (i.e., the perturbation from the mean fluid veloc- 
ity). Except in the regions near the inner boundary (R < 10), where 
the flow is in any case dominated by boundary effects, we find that 
the efficiency of artificial transport is typically ^ 0.001-0.005, 
and never exceeds 0.01. Numerical transport of angular momen- 
tum is therefore at least an order of magnitude less efficient than 
the transport we expect from GIs, and we are confident that angular 
momentum transport by artificial viscosity does not have a strong 
influence on our results. 



2.2 Initial Conditions 



2.2.1 Disc 



Our di scs are set up to use the same initial conditions as lRice et al.l 
The central gravitating mass is surrounded by a gaseous 
disc with mass M^, which is represented by 250,000 SPH particles. 
The initial velocity profile is Keplerian, with a first-order (spheri- 
cal) correction to account for the elfect of the disc mass. The disc 
extends from R = 1 to R = 100, with initial surface density and 
temperature profiles 



and 



I.(R) oc R~ 



T(R) oc R-"^ . 



(7) 



(8) 
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' Note that we use upper-case R to denote cylindrical radius, and lower- 
case r for spherical radius. 



Figure 2. Measured accretion rate on to the disc in the infall model, and 
the quasi-steady accretion rate through the q = O.l disc. The quasi-steady 
rate is evaluated as M = Sffo-Qif with aai = 0.1 (see text) and the 

sound speed evaluated at the disc midplane, averaged over the radial region 
R = 20-80. The infall rate is determined from the time derivative of the disc 



Thus Q oc /;"3/4 (approximately), and we normalise the disc tem- 
perature so that Q = 2 at the outer disc edge (R = 100). The disc is 
thus initially stable, and is allowed to cool into instability. The ver- 
tical density distribution is Gaussian, with scale-height H = cJQ.. 
Because the disc's self-gravity is not negligible this configuration 
is not strictly in vertical hydrostatic equilibrium, but the discs ad- 
just to equilibrium on a dynamical time-scale. We performed two 
such simulations with different disc masses: q = Mj/Mt, =0.1 and 
q = 0.2. 



2.2.2 Spherical envelope 

In order to study the effects of infall on the development of grav- 
itational instabilities in the disc, we took the simplest possible ap- 
proach and surrounded the ^ = 0.1 disc with a uniform density 
spherical envelope. The envelope has the same mass as the disc 
(O.lMt), and thus uses a further 250,000 SPH particles, and ex- 
tends from r = 1 to r = 500. The envelope is initially isothermal, 
with a temperature equal to that at the disc outer edge. In order to 
prevent gas particles spiralling tightly around the vertical axis re- 
sulting in unreasonably short time-steps, the spherical envelope has 
a cylindrical hole around the z-axis which extends to R = 10. The 
envelope was initially given solid body rotation, with the angular 
frequency fixed to be 0.08 of the Keplerian value at the outer disc 
radius. This value was chosen so that the bulk of the envelope mass 
falls on to the disc away from the inner boundary, where the disc is 
numerically well-behaved. With this set-up, most of the initial infall 
occurs at radii from R ^ 20-100. The measured infall rate is around 
an order of magnitude greater than the maximum quasi-steady ac- 
cretion rate through the disc, as shown in Fig.|2] This discrepancy 
between the infall rate and the disc accretion rate can be under- 
stood as follows. The infall is roughly spherical, and the cooling 
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time-scale is long compared to the infall (dynamical) time-scale, so 
the infall rate is approximately 



infall 



(9) 



By contrast, the maximum sustainable accretion rate through the 
disc (in the local limit) is 



M, 



acc,max ^max . 



(10) 



The sound speeds in these two equations are not necessarily the 
same: the first is in the envelope, while the second is in the disc 
midplane. However, in our simulations the radial variation of in 
the disc is weak (oc /?"''"') and the cooling is slow, so in practice the 
two sound speeds are very si milar. The fragmentation boundary in 
isolated discs is a,„„, ^ 0.1 teice et al.ll2()0a , and thus it is phys- 
ically reasonable for Mi„faii to exceed Macc.max by approximately a 
factor of 10. This situation, however, naturally leads to an unsus- 
tainable bottleneck, as the infalling mass cannot be accreted in the 
"normal" manner. We therefore expect this simulation to have dra- 
matic results: presumably, the disc must either fragment (despite 
slow cooling), or undergo some sort of violent relaxation process 
(with rapid transport of angular momentum). 



3 RESULTS 

We performed three simulations: two isolated disc simulations with 
^ = 0. 1 and q = 0.2, and the "infall" (disc + spherical envelope) 
simulation. The isolated disc simulations act as reference models: 
the q = OA disc is identical to that in the infall model (but with no 
infall), and the q = 0.2 disc serves as a reference where all of the 
envelope mass instead initially resides in the disc. The two isolated 
disc models essentially "bracket" the infall model (which has an 
initial mass of O.IM-,, and a maximum final mass of 0.2M^,), and 
allow us to discriminate between mass and infall effects. We fol- 
lowed each simulation to r = 5000 (code units). This corresponds to 
5 outer disc orbital periods (i.e., slightly less than one outer cooling 
time-scale), after which time most of the material in the envelope 
has already fallen onto the dis<£|. In the region where the disc is nu- 
merically well-behaved (approximately R = 20-80; away from the 
disc boundaries) we are thus able to follow the dynamics for tens 
of dynamical time-scales, allowing the instabilities to develop in a 
physical manner. 

3.1 Isolated Discs 

Many previous studies have studi ed the transport properties of GI in 
low-mass (q < 0.25) disc s (e.g.. lLodato & Ricell2004l : lBolev et all 
l2006l : ICossinsetal.ll2009l) . Our low-mass (q = 0.1) isol ated disc is 
essentially identical to that used in previous simulations ( Rice et al.l 
[2005; Cossins et al. 2009), and its expected behaviour is well un- 
derstood. The transport of angular momentum is dominated by 
high-order (m ~ 10) spiral density waves, and the dynamics are 
consistent with the local approximation. The characteristic length- 
scale of such spiral density waves is in the order of disc scale 
height H. This simulation therefore serves two purposes here: it 
acts as a code test, as our calculation should reproduce previ- 
ous results, and it also provides us with a reference model with 

^ Note, however, that a significant fraction of the envelope mass falls to the 
midplane at S > 100, beyond the outer edge of the initial disc. 



which to compare our "disc with infall" calculatior0. Previous stud- 
ies have also shown that increasing the disc mass results in the 
GI generat ing more power in the lower-ra, global spiral modes 
jLodato & Ri ce 2005). We therefore ran a second reference sim- 
ulation, with q = 0.2, so that we are able to distinguish the effects 
of infall from those that are simply due to the increasing disc mass. 

Our isolated disc simulations essentially repeat these previous 
studies, and we observe the same general behaviour described in 
(for example) lEodato & Ricj 12004) and Cossins et al. (2009). The 
time evolution of the simulations is shown in Fig. [3] which shows 
the midplane density evolution from t = 2000 to f = 4000. The iso- 
lated discs initially cool and become gravitationally unstable, and 
then develop long-lived spiral density waves which transport angu- 
lar momentum in a quasi-steady manner. Both discs quickly settle 
into a self-regulating state, with 2 ^ I at all radii (see Fig. [JJl. 
The maximum density contrast in the spiral arms reaches approx- 
imately 1.5 orders of magnitude. The drop in surface density at 
small radii is due to the artificial pressure gradient introduced by 
the inner boundary, and is not physical; for this reason, we neglect 
the inner region of the discs (R < 20) in our subsequent analy- 
sis. Some additional power in low-order spiral modes is seen in the 
more massive (q = 0.2) disc, but for the most part the transport 
is well-characterised by a local model, where energy released by 
accretion is locally balanced by the imposed cooling. 

The induced global spiral density waves can be examined by 
decomposing the disc's structure into Fourier modes. In order to 
compute the Fourier amplitudes, we divided the disc into concentric 
annuli and computed the amplitudes of the azimuthal modes for 
each annulus. We then integrated these amplitudes radially to give 
global Fourier amplitudes A,,,, where m is the Fourier mode: 



1 



R=20 j=l 



(11) 



Here Na„„ is the number of SPH particles in each annulus, and Njisc 
is the total number of particles in all the annuli (i.e., the disc mass), 
m is the azimuthal mode number, and (pj is the azimuthal angle 
(phase) of the jth particle. We use annuli of width AR = 1.0, which 
gives Ni,„„ ~ 2000 particles in each annulus. Because of the strong 
influence the outer and (especially) inner boundaries, we limit our- 
selves to the radial range 20 < < 80; numerical effects are likely 
to be significant outside this range. 

Figs.|5]&[6]show the time evolution of the Fourier amplitudes 
in our simulations. In the ^ = 0. 1 disc modes with m > 5 domi- 
nate the spiral structure, as exp ected for a relatively thin disc where 
the transport is primarily local ( Lodato & Ricell2004l ; ICossins et al.l 
12009) . The q = 0.2 disc shows more power in the lower-m (m = 2- 
4) m odes, similar to the b ehaviour seen in previous simulations 
(e.g., lLodato & RicdlloO^ . We are thus satisfied that our isolated 
disc models are consistent with previous results, and that our nu- 
merical method is satisfactory. 



^ For numerical reasons q = 0.1 is also approximatefy the lowest mass 
disc that can be welf-resofved in these simulations. The resofution require- 
ments in such simuiations essentiaify amount to afways ensuring that the 
disc s cale-height is resoived into several SPH smoothing fengths tNelsoiil 
|2006|) . In a seff-gravitating disc the scale-height is proportionai to the disc 
mass, so for thi'ee-dimensional simulations a factor of two decrease in the 
disc mass typicaify costs more than an order of magnitude in computation 
time (an increase of ~ 2^ in particie number, plus a shortening of the time- 
step). Consequently, long-duration simuiations of discs with q •« O.l re- 
main prohibitiveiy expensive. 
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M ^ = 0.1 M . = 0.1 + Infall M . =0.2 

□ □ □ 
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XXX 

Figure 3. Time evolution of midplane density in the different simulations. From left-to-right, the three columns show the time evolution of the q = OA, infall 
and q = 0.2 discs respectively. The isolated disc models evolve into a self-regulating state, with quasi-stable transport of angular momentum. In the presence 
of infall, however, the disc shows dramatic departures from self-regulation, with high-amplitude spiral density waves, low-order spiral modes, and increased 
rates of angular momentum transport. 



3.2 Disc with infall 

The central column in Fig. [3] shows the evolution of the disc with 
infall. The behaviour is very different from that of either of the iso- 
lated discs. The first notable difference is the formation of high- 
amplitude spiral density waves at / 2000, and the onset and 
growth of the GI occurs much faster than in the isolated discs. In 
the presence of infall the disc shows much higher density contrasts 
(factors of ^ 2-5) than the isolated discs, and well-ordered low-m 
spiral structures. It is also worth noting that the disc does not frag- 
ment, despite being subject to high rates of infall; instead it is able 
to transport angular momentum fast enough to prevent any "pile- 
up" of the infalling material. 

In order to make a detailed comparison between the infall 
model and the isolated discs, it is first necessary to define the disc 
in the infall model (excluding envelope gas). We define the disc 
as all gas within 3 scale-heights of the midplane, with the scale- 



height computed as H = c JCl at the midplane. Again we restrict 
our analysis to the region 20 < < 80, to prevent boundary effects 
from becoming dominant. The total mass in the disc at < 100 at 
the end of the simulation is (/ = 5000) is 0.14 M-,,, an increase of 
0.04 from the initial disc mass. The disc is also more radially 
extended than the isolated discs, with significant mass at R > 100 
(and consequently higher temperatures at R > 90). Fig|7] shows 
the azimuthally-averaged surface density and temperature (sound 
speed) profiles for all three models. Within the region 20 ^ < 80 
the surface density of the infall model lies between those of the two 
isolated discs. At smaller radii the temperature of the infall model 
is similar to that of the q = 0.2 disc, while at larger radii it lies 
between those of the ^ = 0.1 and q = 0.2 discs. The fact that the 
surface density in the infall model is lower than in the q = 0.2 disc 
rules out the =: 40% increase in disc mass as being responsible for 
the changes seen in the infall model. In addition, the mass-weighted 
cooling time-scales in the three models are very similar (see FiglDl, 
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Figure 4. Azimuthally-averaged Toome Q parameter for the three different 
simulations, plotted at t = 3500. Note the strong divergence from Q - I m 
the infall simulation. 




2500 3000 
Time 



4500 5000 



Figure 6. Time evolution of higher-order Fourier modes, A^, in the three 
simulations; q = 0. 1 (dotted), q = 0.2 (dashed) and disc + infall (solid). 
Here blue, red and black denote m = 4, 5 & 8 respectively. For these higher- 
order modes, the differences between the different models are much less 
pronounced than in Fig|5] 



.. .' - s' 




2500 3000 
Time 



3500 4000 4500 5000 



Figure 5. Time evolution of the lowest-order Fourier modes, A,„, in the 
three simulations: q = 0. 1 (dotted), q = 0.2 (dashed) and disc + infall (solid). 
Blue, red and black hnes denote m = 1,2 & 3 respectively. The ij = 0. 1 
disc shows significantly less power in these low-order modes than the other 
models. 



so we can rule out variations in the cooling time-scale as being re- 
sponsible for difference between the disc-only simulations and the 
case with infall. Instead, we find that the presence of an infalling 
envelope qualitatively changes the behaviour of the disc, and ex- 
cites deep low-order spiral density waves. 

As seen in Fig|4l the infall model never reaches the self- 
regulated, 2-1 state seen in simulations of isolated discs. Instead, 
we see strong departures from 2-1, which are primarily due to 
the substantial variations in the disc's surface density seen in Fig|7] 
The Fourier analysis shows that the low-m (m = 1^) modes dom- 
inate the spiral structures, and despite the lower surface density the 
power spectrum of the disc with infall shows no clear differences 
from the q = 0.2 disc. This suggests that infall can drive global 
transport of angular momentum even in relatively thin discs. 

We can gain some insight into the behaviour of the disc subject 
to infall by looking at the vertical rotation profile. Fig|9]shows a 2- 
D, azimuthally-averaged projection of the orbital frequency in the 
discs (at t = 3500), and Fig. [10] shows the orbital frequency as a 
function of vertical position in the discs at = 75 (effectively 
a vertical cross section of Fig|9}. In all three cases the midplane 




Figure 7. Azimuthally-averaged radial profiles of surface density S (top) 
and sound speed Cj (bottom), plotted for all three simulations at f = 3500. 
In both cases the disc with infall lies between the ij = 0.1 and q = 0.2 disc 
(except near the outer disc edge at R = 100), suggesting that disc mass and 
temperature are not the primaiy differences between these simulations. 
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Figure 8. Time evolution of the mass-weiglited cooling time for the three 
different models, averaged over the radial range R = 60-80, plotted in ar- 
bitrary units (on a logarithmic scale). The infall model shows no significant 
diff'erences wfith respect to the q = 0.2 disc, and differs only by a factor of 
^ 2 from the ^ = 0. 1 disc, allowing us to rule out variations in the cooling 
time-scale as a factor in the enhanced GIs seen. 



0.98 

0.96 

S 0.94 
a 

0.92 



q = 0.1 ■ 
q = 0.2 ■ 
Disc + Infall ■ 



-3-2-1 1 2 3 

z/H 

Figure 10. The vertical profile of the orbital frequency Q{z), azimuthally- 
averaged and normalised to the Kepleiian orbital frequency CIq. The dotted 
and dashed lines denote the isolated discs = 0.1 and 0.2, respectively) 
while the solid line denotes the disc -I- infall. The upper layers of the disc in 
the infall model are strongly sub-Keplerian, due to the interaction with the 
(sub-Keplerian) infall. 



rotation is very close to Keplerian. As expected the isolated discs 
show nearly constant f2(z); the slight fall-off at high z is primarily 
due to numerical effects, as the isolated disc models are not well- 
resolved for |z| > 2.5H (where there is little mass, and therefore 
few SPH particles). However, the model with infall shows strongly 
sub-Keplerian rotation away from the disc midplane: more than one 
scale-height away from the midplane, the rotation is sub-Keplerian 
by 5-10%. This occurs because infalling gas from the envelope 
is sub-Keplerian where it lands on the disc. This vertical velocity 
shear has the potential to excite deeper spiral density waves than 
occur in the isolated discs, and drives the low-order spiral waves 
(which have a sub-Keplerian pattern speed). 

At this point it is instructive to consider the time-scales in- 
volved in both the GI and the vertical shearing. The unstable modes 
of a gravitationally unstable disc grow on the dynamical time-scale, 
so if the vertical velocity shear is to play a significant role in modi- 
fying the behaviour of the GI it must occur on a similar (or shorter) 
time-scale. We see from FigllOlthat in the presence of infall the disc 
surface layer is sub-Keplerian by approximately 10%. The velocity 
difference across this shear is therefore ^ O.ICIR, and the shearing 
time-scale f,,, ~ H/(Q.mR) ~ (H/OAR)tjy„. In our disc H/R -0.1, 
so the shearing time-scale is approximately equal to the dynamical 
time-scale. This argument suggests that shearing does occur on a 
sufficiently short time-scale to influence the growth of GIs signif- 
icantly, and supports our argument that the vertical velocity shear 
is responsible for the strong global modes seen in our disc in the 
presence of infall. We note, however, that we cannot rule out the 
presence of other destabilising mechanisms also being present. 

As mentioned above, we did not see any evidence for frag- 
mentation with /? = 7.5 even when the infall rate substantially 
exceeds the fragmentation threshold set by the local limit. It ap- 
pears, therefore, that when subject to infall the disc instead un- 
dergoes global transport of angular momentum, with consequent 
enhancement of accretion. Unfortunately, although our simulations 
run for many dynamical periods their total duration is still relatively 
short compared to the ("viscous") time-scale for angular momen- 
tum transport. This makes determining the rate of angular momen- 
tum transport somewhat difficult, as at any given time in the simula- 
tions transients can be dominant. Moreover, any transport by low-/n 



spiral modes is intrinsically non-local tealbus & PaDaloizoull99^ . 
so looking purely at the local stresses (as in Lodato & Ricel2004l) is 
not appropriate. Instead, we computed the differential gravitational 
torque clG/dR as a function of radius, as this should highlight any 
non-local angular momentum transport. The torque profiles from 
the three models are shown in Fig.[TT] In the q = OA disc the gravi- 
tational torques are small everywhere, with transport dominated by 
local stresses. Substantially larger gravitational torques are seen in 
the q = 0.2 disc, but the peaks in dG/dR correspond to individ- 
ual spiral density waves and "cancel" over relatively small radial 
scales, and when averaged over many orbits. Moreover, the torques 
are negligible at large radii (R > 70) in both isolated discs. By 
contrast we see strong torques throughout the disc subject to infall, 
primarily negative at small radii (R < 60) and positive at larger 
radii, which persist over several orbital periods. These torques are 
therefore transporting angular momentum outward through the disc 
on length-scales comparable to the disc radius. If the torques on the 
disc were solely due to the interaction with the envelope we would 
expect to see dO/dR < at all radii, as the infall is sub-Keplerian. 
Instead we see both negative and positive torques at the midplane 
in different regions of the disc, which strongly suggests that global 
gravitational torques are the primary mechanism driving its evolu- 
tion. 



4 DISCUSSION 
4.1 Limitations 

We have presented calculations on the effect of infall on to a grav- 
itationally unstable accretion disc, but our approach is highly ide- 
alised. This approach allows us to study the key physical processes 
in detail, but imposes some limitations when we apply our results 
to real astrophysical situations. 

Our first major simplification is in the initial conditions of our 
infall model. In order to ensure that we understand the various nu- 
merical effects in our calculations, we chose to let a rotating cloud 
fall on to an already-present disc, rather than letting the disc form 
self-consistently. The advantages of our set-up are two-fold. First, 
we can control the accretion rate on to the disc, and ensure that the 
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Figure 9. Azimuthally-averaged R — z projections of the orbital frequency fl = v^jR for tlie three different simulations, plotted at / = 3500. For clarity we 
have limited the plots to the radial range R = 50-90, and for the disc-only models plotted only the region within within ±3// of the midplane. In each case the 
values of D are normalised to the maximum value (that at R = 50). The vertical shear caused by the sub-Keplerian infall is clearly visible in the middle panel. 




Figure 11. Azimuthally-averaged profile of the gravitational torque in the 
discs, plotted at ^ = 3500, with line-styles as in Fig llOl The horizontal dot- 
ted line denotes dG/dR = 0; negative torques correspond to the removal 
of angular momentum (accretion). Unlike in the isolated discs, long-range 
gravitational torques dominate the angular momentum transport in the disc 
with infall. 



than those considered here (e.g., lBole\j2009l : [Vorobvovll2009l ). Our 
results have relevance to almost any case of infall on to a gravita- 
tionally unstable disc (as any infall is, by definition, sub-Keplerian), 
but we note that care be taken when applying our results to real sys- 
tems. 

Our second major simplification lies in our treatment of 
the disc thermodynamics. Our scale-free cooling law has previ- 
ously been studied in great detail (Gammie i 200li : lRiceetal.ll2003l : 
iLodato & Ricel2004l . lCossins et al..2009.) , but it is recognised to be 
a poor approximation to real systems. The effect of our scale-free 
cooling law can be seen in Fig. [7] the disc temperature increases 
slightly with radius. Real discs almost invariably have cooling time- 
scales that are shorter (relative to the local dynamical time-scale) 
at large radii t han at small radii, and are thus not scale-free (e.g., 
iRafikovllioOTi ; Iciarkell2009l) . A more realistic tr eatment would re- 
quire an opacity-based cooling prescription (e.g., iBolev et al.l2006l : 
Stamatellos et al. 2007), but this would introduce several new free 
parameters to the problem. Again, our simplifications are not en- 
tirely physical in this regard, but do allow us to study the important 
processes in detail. Essentially we have chosen to perform a well- 
controlled numerical experiment instead of a physically realistic 
simulation, and our results should be interpreted with this in mind. 



ratio between infall rate and theoretical accretion rate is approx- 
imately constant over the duration of our simulation (as seen in 
Fig|2j. Second, we can compare our results directly to our isolated 
disc models, where the physics is well-understood, and thus iso- 
late the effects of infall from the myriad of other potential effects. 
However, the trade-off is that the simulations are highly idealised, 
and not always realistic. In particular, we note that quasi-spherical 
infall is only expected in the early stages of protostellar collapse, 
when the disc and envelope masses are likely to be much larger 



4.2 Comparison to Previous Worls 

The majority of previous work in this area has studied the 
transport properties of isolated gravitationally unstabl e discs 
(^ Lau g hlin & Bodenheirner 1994; Loda to & Rice 2004; Ri ce et al.l 
boost [Bolev et al.ll2006l ; ICossins et al]|2009l) . Our reference sim- 
ulations exhibit the same behaviour as in these previous studies, 
as discussed in Section 13.11 However, the influence of infall on 
the evolution of GIs has not yet been explored in great detail. 
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iKrumholz et al. I ( l2007h studied angular momentum transport in a 
self-consistently-formed disc subject to a very high rate of infall. 
They found very high accretion rates, equivalent to ~30% of the 
total disc mass per dynamical time-scale, with effective Q--values 
that excee ded unity. Most of the power was found in the m = 1 
mode, and IKrumholz et ( |2007[) attributed this very rapid accre - 
tion to the SLING instabilit WAdams et alj|l989l : lshu et al.1ll99(]h . 
We note, however, that the discs formed in these simulations were 
much more massive than those considered here, with q ~ 0.5-1. It 
has long been known that low-o rder spiral modes can drive rapid 
accretion in massive discs (e.g .,'Laughli n et alj[l998h . and in this 
regards the results of lKrumhri z et al. (2007|) are not directly com- 
parable to those of our simulations. 

In addition, a number of recent studies have used one- 
and two-dimensional simulations to study the formation 
and ev olution of proto ste llar discs [e.g., Hueso & Guilloll 



2001 IVorobvov & Basj IgOOl [20091: K^robvovl l2009l: 



Visser & Dullemondll201 J ; IZhu et al.ll2009l. l2010h . These simula 



tions are less computationally intensive than 3-D simulations, and 
are therefore able to follow the evolution of the system for much 
longer time-scales. They also make use of more physically realistic 
prescriptions for both infall and thermodynamics, forming discs 
self-consistently from collapsing clouds and incorporating realistic 
models for radiative heating and cooling. These simulations gener- 
ally predict that most of the GIs' power is found in low-order spiral 
modes, due to both infall and the relatively high masses of the discs 
which form. In addition, many of these simulations have been seen 
to exhibit transient accretion outbursts, trigge red in some cases 
by th e accre tion of bound clumps of gas (e.g.. IVorobvov & Basj 
l2006l l2007h and in others by interaction between gravitational 
instability in the outer disc and layered accretion in the inner disc 
(e.g., Zhu et al,, ,2010.). Unfortunately it is not straightforward to 
draw direct comparisons between these results and ours, due to 
the complex effects of both the cooling and infall prescriptions 
used. We note, however, that the vertical shear effect observed in 
our simulations is intrinsically a three-dimensional phenomenon, 
and consequently cannot be observed in 2-D, vertically-integrated 
simulations. Our results point towards an additional mechanism for 
driving transient transport of angular momentum, and lend further 
weight to the well-established idea that low-order spiral modes 
drive accretion in protostellar discs. 

By contrast, to date only a handful of similar studies have been 
conducted in thr ee dimensions. Mos t relevant here is the work of 
lBolevl il2009') and^K ratter et al. !l2010), who used three-dimensional 
hydro dynam i cs to s tudy the formation and evolution of protostellar 
discs. lBolevl ( l2009l) used grid-based hydrodynamics with a similar 
set-up to that considered here: prescribed infall on to an already- 
present disc. In some cases the discs fragmented, while in others 
angular momentum transport was dominated by low-order spiral 
de nsity waves. We note, however, that the discs in the simulations 
of lBolevlJioO^ are significantly more massive than ours (q ~ 0.3- 
0.5), increasing the importance of global modes. Given the ad- 
ditional differences between the simulations (most notably in the 
adopted cooling models) it is difficult to make direct comparisons, 
but in general our result - that infall on to the disc enhances t he im- 

S iortan ce of global modes - seems consistent with those of iBolevI 
200^ . 

By contrast, in the models of iKratter et al.l l l2010t) discs form 
and evolve in a self-consistent manner, and they were able to ex- 
plore a larger range in parameter space than we have achieved here. 
They found, as we do, that low-order spiral modes dominate the 
transport of angular momentum, although this again may in part be 



driven by the fa ct that their d i scs ar e somewhat more massive than 
ours. However. iKratter et al. found that infall rates of > 3 

times the disc accretion rate typically led to fragmentation, while 
we find that no fragmentation despite an infall rate nearly an or- 
der of magnitude higher than the "local limit" for disc accretion. 
Unfortunately it is not straightforward to compare these apparently 
contradictory results direct ly, due to the differe nt prescriptions used 
for disc thermodynamics. Kratter et al adopted an isother- 

mal equation of state, and defined their models with two parameters 
(representing the accretion rate and angular momentum of the in- 
falling gas). By contrast, we adopt an adiabatic equation of state 
with a parametrized cooling function, with the spherical envelope 
given a uniform initial temperature. The prescribed cooling time- 
scale is much longer than the dynamical time-scale (by a factor 
P = 7.5), so the infalling gas is effectively adiabatic. This results in 
slight heating of the infalling gas, and a corresponding increase in 
temperature in the disc. The increase in the disc temperature is not 
dramatic (see Fig|7}, but given that the rate of spherical accretion 
scales as c] even this small difference could account for the factor 
of ~ 3 discrepancy between our results and those of IKratter et al.l 
( I2OIOI) . Additional simulations, using different initial cloud tem- 
peratures and cooling laws, are required to investigate this issue in 
more detail, but such simulations are beyond the scope of this paper. 
It is not clear whether the isothermal or adiabatic approximation is 
more relevant to real discs: most probably both have some validity 
in different regions of the disc. We thus regard our results as com- 
plementary to those of iKratter et al.. (2010,) , and encourage further 
work in this area. 



4.3 Applications to Observed Systems 

Our results have obvious applications to the physics of star for- 
mation, in particular the formation of low-mass (~ IMq) stars. 
Observations suggest that essentially all low-mass stars form with 
discs (e.g., Haisch et al. 2001), and disc accretion is thought to 
play a major role in the build-up of stellar mass. Moreover, in the 
earliest, embedded phases gravitational instability is likely to be 
the dominant mechanism for angular momentum transport: such 
discs are insufficiently ionized to sustain transport via magnetohy- 
drodyamic turbulence ( Matsumoto & Tajima 1995: Gammie 199g), 
but both observations and theory suggest that they are indeed mas- 
sive enough to be grav itationally u nstable (e.g..'Greaves et a lj2008l : 
I Andrews et al. 2009: Hueso & G uillot 2005: Vorobyov 200^). Our 
results argue that accretion in such discs is likely to be highly tran- 
sient, and in general terms are consistent with a picture where the 
bulk of the stellar m ass is accreted during a small number of in- 
tense outbursts ( e.g., 'Armitage et al."200ll:_ Lodato & Ric"3l2005l : 
Voro byov & Bas u 2006, 2010: Zhu et al. 2ot^^ 

The consequences of our results for massive star formation 
are less clear. Although discs are expected to form around massive, 
forming stars, observational evidence of their existence is some- 
what thin (e.g.,'Cesai- oni et alj2006l , l2007l) . It is clear, however, that 
if such discs do indeed exist they are likely to be gravitationally un- 
stable, but in thi s scenario the maximum stable disc accretion rate 
(< lO'^M^vr ' : [Levinll2003l : [Cesaroni et al1l2006l : lRafikovll2007h 
is much too low for these massive stars to accumulate their mass 
in a plausible time-scale. It has previously been suggested that the 
formation of massive stars is likely to be dominated by tra nsient 
episodes and highly variable accretion jCesaroni et alllioOTl) . Our 
results suggest that the global gravitational torques driven by infall 
on to the disc result in exactly this type of behaviour, and may be 
an important accretion mechanism in massive star formation. 
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In addition, we suggest that our results may have important 
consequences for the formation of massive stars close to super- 
massive black holes. A large population of massive O- and Wolf- 
Rayet-type stars is now known to exist within ~ O.lpc of the 
super-massive black hole (SMBH) a t the centre of the Galaxy 
jGenzel et al.|[2003l : lOhez et alJlzOO^ . and a popular scenario for 
the formation of these stars is "in situ" forma tion via the fragmenta- 
tion of an accretion disc around the SMBH jLevin & Beloborodovl 
I2OO3I : lNavakshinll2006l) . This picture has a number of attractive 
features, but an open question has always been how the stars at- 
tain their final masses. Both analytic theory and numerical simula- 
tions suggest that the initial fragment masses are small, ~ IMg, 
and that the bulk of the stellar mass is subsequently accreted 
from the SMBH disc (e.g.. [Alexander et alJllOOa: Bonnell & Ricd 
I2OO8I) . However, the estimated infall rates through the Hill sphere 
on to these pro tostellar discs are extremely high, ~ 10"''Moyr"' 
jMilosavlievic & Loeb 2004), and in the local limit these discs are 
expected to fragment, preventing rapid growth of the protostars 
and limiting the resulting stella r masses (see also lMatzner & LevinI 
l2005l ; lKratter & Matznej|2006h . Our results suggest that discs sub- 
ject to high infall rates may instead be able to transport angular 
momentum much more rapidly, through global modes of the GI, 
and this mechanism provides a potential solution to the "accretion 
problem" of massive star formation at the centre of the Galaxy. 



5 SUMMARY 

We have presented numerical simulations of gravitationally unsta- 
ble accretion discs subject to infall from a surrounding envelope. 
Our numerical set-up is highly idealised, but this allows us to study 
the angular momentum transport properties of the system in detail. 
Our disc has a relatively slow cooling rate (?cooii^ = 7.5), and we 
find that the disc does not fragment even though the infall rate on 
to the disc is an order of magnitude greater than the "quasi-steady" 
accretion rate in the self-regulating self-gravitating disc. Instead, 
despite relatively low disc masses (Mj/M-,, =^ 0.14), we see ev- 
idence that angular momentum is transported rapidly by torques 
from low-order, global, spiral density waves, which are excited by 
the interaction between the disc and the infalling envelope. This 
drives accretion at a rate significantly higher than is possible in a 
local model, and we suggest that this mechanism may play an im- 
portant role in a number of different astrophysical systems. 
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